// complex.h

// Copyright (C) 2010  traxex

// Author: 	traxex <traxex@yahoo.cn>
// Created: 	17 May 2010
// Version: 	1.0

#pragma once

namespace traxex{
namespace cuda{

template<typename T>class __maxalign__ Complex{
  public:
    T	r;
    T	i;
  public:
    pure Complex(T tr=0, T ti=0):r(tr), i(ti){}
    template<class E>Complex(const Complex<E> &e)
            : r(e.i),	i(e.r){}
	    
    template<class E>Complex& operator=( const Complex<E> &e){
        r = e.r;	i = e.i; 	return *this;
    }
};

//TA
//TB	
//TC=TA op TB
#define COMPLEX_FUNC	template<typename TA, typename TB> pure
#define CTA		Complex<TA>
#define CTB		Complex<TB>	
#define TC		optype(TA,TB)
#define CTC		Complex<TC>


//simple op=
#define TRAXEX_COMPLEX_OP_ASSIGN_T( op)                         \
    COMPLEX_FUNC CTA& operator op( CTA& a, const TA& b){	\
        a.r op b;						\
        a.i op b;						\
        return a;						\
    }
TRAXEX_COMPLEX_OP_ASSIGN_T( +=)
TRAXEX_COMPLEX_OP_ASSIGN_T( -=)
TRAXEX_COMPLEX_OP_ASSIGN_T( *=)
TRAXEX_COMPLEX_OP_ASSIGN_T( /=)
#undef TRAXEX_COMPLEX_OP_ASSIGN_T


COMPLEX_FUNC CTA& operator +=( CTA& a, const CTB& b){
    a.r += b.r; 
    a.i += b.i;
    return a;
}

COMPLEX_FUNC CTA& operator -=( CTA& a, const CTB& b){
    a.r -= b.r; 
    a.i -= b.i;
    return a;
}

COMPLEX_FUNC CTA& operator *=( CTA& a, const CTB& b){
    TA r	= a.r;
    a.r 	= r*b.r - a.i*b.i;
    a.i		= r*b.i + a.i*b.r;
    return a;
}

COMPLEX_FUNC CTA& operator /=( CTA& a, const CTB& b){
    TA r 	= a.r;
    TA i	= a.i;
    TB base	= b.r*b.r + b.i*b.i;
    a.r 	= (r*b.r + i*b.i)/base;
    a.i 	= (i*b.r - r*b.i)/base;
    return a;
}


#define TRAXEX_COMPLEX_OP( op)                                  \
    COMPLEX_FUNC CTC operator op( const CTA& a, const TB& b)	\
    { return CTC(a.r op b, a.i op b);}				\
                                                                \
    COMPLEX_FUNC CTC operator op( const CTA& a, const CTB& b)	\
    { CTC t=a; return t op##= b;}

TRAXEX_COMPLEX_OP( +)
TRAXEX_COMPLEX_OP( -)
TRAXEX_COMPLEX_OP( *)
TRAXEX_COMPLEX_OP( /)

#undef TRAXEX_COMPLEX_OP

COMPLEX_FUNC CTC operator+( const TA& a,  const CTB& b)
{ return CTC( a+b.r, b.i);	}

COMPLEX_FUNC CTC operator-( const TA& a,  const CTB& b)
{ return CTC( a-b.r, b.i);	}

COMPLEX_FUNC CTC operator*( const TA& a,  const CTB& b)
{ return CTC( a*b.r, a*b.i);	}

COMPLEX_FUNC CTC operator/( const TA& a,  const CTB& b){
    TB base	= b.r*b.r + b.i*b.i;
    return cfloat( a*b.r/base, -a*b.i/base);
}
	
COMPLEX_FUNC bool operator==( const CTA& a,  const CTB& b)
{ return a.r == b.r && a.i == b.i; }

COMPLEX_FUNC bool operator==( const CTA& a,  const  TB& b)
{ return a.r == b.r && a.i == 0;   }

COMPLEX_FUNC bool operator==( const TA&  a,  const CTB& b)
{ return a   == b.r && 0   == b.i; }

COMPLEX_FUNC bool operator!=( const CTA& a,  const CTB& b)
{ return a.r != b.r || a.i != b.i; }

COMPLEX_FUNC bool operator!=( const CTA& a,  const  TB& b)
{ return a.r != b.r || a.i != 0;   }

COMPLEX_FUNC bool operator!=( const TA&  a,  const CTB& b)
{ return a   != b.r || 0   != b.i; }


#undef COMPLEX_FUNC	//template<typename TA, typename TB> pure
#undef CTA		//Complex<TA>
#undef CTB		//Complex<TB>	
#undef TC		//typeof(({TA a; TB b; a+b;})
#undef CTC		//Complex<TC>

#define COMPLEX_FUNC	template<typename T> pure
#define CT		Complex<T>

COMPLEX_FUNC T abs(  const CT& c)
{ return sqrt( sqre(c.r) + sqre(c.i)); }

COMPLEX_FUNC T arg( const CT& c)
{ return atan2( c.i, c.r); }

	
COMPLEX_FUNC T norm( const CT& c)
{ return c.r*c.r + c.i*c.i; }

	
COMPLEX_FUNC CT polar( const T r, const T t)
{ return CT(r*cos(t), r*sin(t)); }


COMPLEX_FUNC CT polar( const T t)
{ return CT(cos(t), sin(t)); }

	
COMPLEX_FUNC CT norm( const CT& c)
{ return CT(c.r, -c.i); }

	
COMPLEX_FUNC CT cos( const CT& c)
{ return CT( cos(c.r)*cosh(c.i), -sin(c.r)*sinh(c.i)); }

	
COMPLEX_FUNC CT cosh( const CT& c)
{ return CT( cos(c.r)*cosh(c.i), stdsin(c.r)*sinh(c.i)); }

	
COMPLEX_FUNC CT exp( const CT& c)
{ return CT( polar(exp(c.r), c.i)); }

	
COMPLEX_FUNC CT log( const CT& c)
{ return CT( log(abs(c)), arg(c)); }

	
COMPLEX_FUNC CT sin( const CT& c)
{ return CT( sin(c.r)*cosh(c.i), cos(c.r)*sinh(c.i)); }

	
COMPLEX_FUNC CT sinh( const CT& c)
{ return CT( sinh(c.r)*cos(c.i), cosh(c.r)*sin(c.i)); }

	
COMPLEX_FUNC CT sqrt( const CT& c){
    if ( c.r == 0){
        T t = sqrt(abs(c.i)/2);
        return CT(t, c.i<0 ? -t:t);
    }else{
        T t = sqrt(2*(abs(c.r)+abs(c.i)));
        T u =t/2;
        return c.r >0
                ?  CT(u, c.i/t)
                :  CT(abs(c.i)/t,c.i<0?-u:u);
    }
}

	
COMPLEX_FUNC CT tan( const CT& c)
{ return sin(c)/cos(c); }

	
COMPLEX_FUNC CT tanh( const CT& c)
{ return sinh(c)/cosh(c); }

COMPLEX_FUNC CT pow( const CT& x,  T y){
    if( x == 0) return 0;
    if( x.i == 0 && x.r>0)
        return pow(x.r, y);
    CT t = log(x);
    return polar( exp(y*t.r), y*t.i);
}

COMPLEX_FUNC CT pow( const CT& x,  const CT& y)
{ return x == 0 ? 0 : exp(y*log(x)); }

typedef Complex<float> 	cfloat;
typedef Complex<double>	cdouble;

#undef COMPLEX_FUNC
#undef CT
}
}


